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^yj . Abstract: Reaction diffusion systems with Turing instability and mass conservation are 

• I studied. In such systems, abrupt decays of stripes follow quasi-stationary states in sequence. 
H , 

At steady state, the distance between stripes is much longer than that estimated by linear 
stability analysis at a homogeneous state given by alternative stability conditions. We show 
^ ' that there exist systems in which a one-stripe pattern is solely steady state for an arbitrary 

oo ■ 

size of the systems. The applicability to cell biology is discussed. 
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Pattern formation, the emergence of a spatial structure from an initially uniform state, 
has been often studied on the framework of reaction diffusion systems (RDS). It is exten- 
sively applied tophysical, chemical, and biological systems to explain their specific spatial 
structures 0,0, 3|. Turing instability is the most prominent mechanism, forming spatially 



periodic stripes 



^. The intrinsic distance between stripes is, in principle, estimated by the 



linear stability analysis at a homogeneous state 



2,|3|,y]. However, this estimation could 



be invalid when applied far from a uniform state. For example, a second bifurcation can 
arise which would indicate the collapse of a simple periodic structure Pi . In such situations, 
the transient dynamics of pattern formation would be difficult to predict. In general, RDS 
shows various dynamics even when steady states are reached j|6|. So far, few studies have 
discussed transient dynamics using the computational analysis of the famous Gray-Scott 
model [6,Q] and by reduced dynamics on the slow manifold jsl O. 11^. 

In this Letter, we study a class of RDS in the context of the above aspects. We consider 
RDS showing Turing instability, in which no production and no degradation of substances 
occur As we will see, the following properties are observed in common; (i) The transient 
dynamics is a sequential transition among quasi-steady states, with a decrease in the number 
of stripes, (ii) The distance between resultant stripes cannot be estimated from the linear 
analysis at uniform state. In particular, there are systems in which a one-stripe pattern is 
solely stable state regardless of the system size. 

Consider a diffusible chemical component with two internal states, U and V. Diffusion 
coefficients are and D^, respectively, for which we can set Du < without loss of 
generality. The transition rates between U and V are regulated by each other. We studied 
a 1-dimensional system with size L(0 < x < L) under periodic boundary conditions unless 
otherwise stated. Concentrations in U and V at position x and at time t are represented by 
u{x,t) and v{x,t) respectively, and obey the following equations. 

dtu = Dudlu - f{u, v) (1) 
dtV = D,dlv + f{u,v) (2) 
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Obviously, the total quantity of the substances (total mass) is conserved. 
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L 



s = Y i'^ ~^ v)dx. (3) 
L Jo 

s is the average concentration of the substance, which is determined by the initial condition 

u{x, 0) and v{x, 0). 

Below, all numerical simulations were performed with f{u,v) = au/{b + m^) — v where 
s = 2.0, a = 1.0, b = 0.1, Du = 0.02, and D„ = 1.0. We observed qualitatively the same 
phenomena in several mass-conserved models 12|. 

Uniform state w* = {u*,v*) is derived from following the conditions; u* + v* = s and 
f{u*,v*) = (stable fixed point in kinetic equation). Let /* (/*) be partial derivatives of / 
with regard to u (v) at w*. If the following relations are satisfied, uniform state w* loses its 
stabihty and the pattern starts to rise. 

/: < /: < 0, (4) 
dj: - dj: > 0. (5) 

All the waves (e*'^^) with wave number k between < < {D^f* — D^f*) / D^Dy are 
unstable. At the beginning of the dynamics, the wave with the largest instability grows 
as in usual Turing systems (see A in Fig. ^ and the line segment representing the most 
unstable wavelength £„,)• 

In a mass-conserved system, characteristic transient processes are observed. After the 
growth of a number of stripes (A in Fig. some stripes stop growing and begin to decay 
(B). With the decay of a stripe, neighboring stripes grow due to mass conservation. The 
distance between neighboring stripes becomes larger (B-C). If the distance is large enough 
the state appears to reach a steady state (C, quasi-steady state). However, one (or more) 
stripe(s) collapses abruptly with the concomitant growth of adjacent stripes (D). As the 
process continues, the number of stripes decreases and the intervals between the abrupt 
transitions gets longer (notice the log-scale representation in Fig. P). In Fig. ^ the system 
finally reaches a one-peak state. The wavelength is much larger than im- Similar processes 
were observed in many mass-conserved systems. 



To understand the observed transient processes, consider the stationary patterns of the 
system with size L. A stationary pattern 'Wq{x) = {uq{x) , Vq{x)) is given by the solution of 
Eq. pi2|) with left hand sides replaced by 0. In a mass-conserved system, there is a family 
of stationary solutions parameterized by s, represented by wo{x; s) explicitly. A function 
h{x) and a value P exist, such that Uq = h{x)/Du and Vq{x) = {—h{x) + P)/Dy, and satisfy 
the following equations. 



Notice DuUo{x)+DyVo{x) is independent of x. P is related to s by P = DyS — {D^—Du)h/D. 



Let us represent the linear operator at a stationary state wohj C. There are two eigen 
functions belonging to the zero eigen value (0-eigen functions); d^w^ = {d^UQ^dxVo) and 
dsWo = {dsUo,dsVo). The former function is derived from the fact that the arbitrary 
translation of stationary state, wqIx + A) is also stationary, while the latter is from the 
conservation property of the mass. Conjugate operator C* is defined as the transposed 
matrix of C One of the 0-eigen functions of this operator is (f){x) = (1, 1). The inner 
product between / = {f^,f^) and g = {g^,g^) (each defined on < x < L) is to be 
{fi9)L = x Joif^g^ + f^9^)dx. Then, is thought to be the conjugate vector of dsWo 
because (0, dgWo) l = 1 and (0, dxWo) l = 0. 

Now, to evaluate the stability of a stationary pattern, consider the following situation. 
Take the one-stripe stationary state Wq in the system with length, which takes the min- 
imum uq aX X = 0(= ^) and the maximum at x = Then copy the exact same state 
on ^ < X < L, and name the system on < x < L as the non-perturbed system (NPS). 
Left and right halves are independent of each other. Next, the boundary condition in NPS 
is changed at x = 0(= ^) and L) into the usual periodic boundary condition of the 
system on < a; < L. We refer to this modified system, which is the one we are interested 
in, as the perturbed system (PS). We represent the state constructed as above by Wq © Wq, 
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P = DuUo{x) + D^vq{x) 

d?h(x) _ f ( _h_ -h+p \ 
<ia;2 J \Du' ) 



(6) 



(7) 



in which h = h{x)dx is the average of h{x). 



where the left (right) hand side of © represents the function on < x < ^ < x < L). 
This state is obviously a stationary solution in both NPS and PS. 

Linear operators at the state are given by Cq for NPS and C for PS. Because NPS is 
simply the juxtaposition of identical systems, ip'l = d^Wo © dxtvo and ip2 = dgtvo © dgWo are 
0-eigen functions for Cq. They are also 0-eigen functions in PS; Cqi/j^ = Cip^ = (z = 1, 2). 



'unctions of Cq is Tp^^ = dsWo © {—dgWo) but this is not 0-eigen functions 



Another 0-eigen 

of C anymore [1^. However, when the amplitudes of \dsUo\ and \dsVQ\ are small at x = 
and ^, the discrepancy between NPS and PS is small and we can expect a function ip and a 
value A which are close to and 0, respectively, and that they satisfy following relation. 

Cif = Xip (8) 

If A is positive, then the stationary state Wq © is unstable and small fluctuations grow 
as ~ e^^ip- Because ip is similar to ip^ = dsWo © {—dsWo), the corresponding dynamics 
appears as the decay of a stripe and the growth of the other, as is observed in the numerical 
simulations. Note the conjugate function of V's is 0l = © (— 0)- 

To check the validity of the above considerations, we numerically measured some related 
quantities. At first, we simulated the equations in the -^-length system and obtained a steady 
one-stripe state. Then, we extended the system size twice and copied the steady state to the 
region k < x < L. Next, perturbations were added keeping the total mass conserved and 



IS 



observed the resulting dynamics. |As(t)| = K^l,"?^)/,! = j {u + v)dx — Jl{u + v)dx 
plotted in Fig. |2l^a) which shows the exponential growth of the perturbation. Then, we 
compared ip^ with the growing part of {u,v). In Fig. I2tb), Aw{x) = {Au{x), Av{x)) is 
shown, where Aw(x) = ^2) — w{x,ti) is the difference of w between two growing time 
points ti and ^2- Aw{x) is similar to ip^ which validates the above considerations. 

The expected '0 is a continuous and smooth function on < x < L and odd around 
X = ^. Thus, it is enough to consider a nontrivial solution of Eq. (jH)) on < x < with 
boundary condition '0(0) = 0'(-|) = 0. We can limit our arguments on < x < in the 
following discussion. We properly redefine C under this limitation. 
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To obtain ip and A in Eq. (jH)), ^ = dgWo + ff is defined. Here, ff = {rj^, r]^) is orthogonal 
to dsWo, i.e. (0,^L/2 = 0. In the first order of approximation, fj satisfies the relation 
Cf] = XdsWQ. The terms m{x) and Q^{x) are introduced in ri^{x) = m{x)/Du and Vj^ = 
(— m(x) + Q^{x))/ Dy. The following equations are then obeyed. 

^ = X{dsUo + d,Vo) (9) 



& - - ft) = ft^H^) + A5.no (10) 



From Eq. © and the boundary conditions, Q{x) = DudgUo + DydgVo + Q^{x) is given by 
Q{x) = \Q{x) where 

f-x nx' T 

Q{x) = / dx' dx"idsUoix") + dsVoix")) - -x (11) 
Jo Jo 4 

By substituting Q^{x) into Eq. (fTUI). m(x) = \A{x) — S^P is obtained, where A{x) is the 

solution of the following equation with the boundary condition A{0) = A{^) = 0. 

Due to the orthogonality {(p,ff)L/2 = 0, A is given by 

X=(^^^—£^A + ^q] \ (13) 



^ -DuD-^ Dy ^ 

where A = J^^ A{x)dx and Q = J^^ Q{x)dx. 

We calculated Eq. (jl3|) numerically. Observed growth rates and estimated values from 
Eq. (fT^ are plotted in Fig. 121 (a) against half of the system size -j, which is the distance 
between two stripes. The two plots are in good agreement with each other and support the 
validity of the arguments presented. 

From Q{x) = AQ(x), we can obtain A = -f^- In NFS, Q(x) is given by 9,P©(-(9,P) 
and independent of x, while in PS Q{x) is connected with at x = and = -I by the 
modification of boundary conditions. Thus, the sign ofM2) is the same as dgP and the sign 
of A is the opposite of dgP. This leads to the condition of two stripe instability in L-length 
system as 

dsP < 0. (14) 
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The above arguments with the two stripes situation are extensible to an identical A^- 
stripe pattern, where each stripe has width. Consider a set of independent functions 
^0 — ©jLi e*^''-'<9st?o = 1)2,--- ,A^), where is defined similarly to © and Wq is 

redefined by the one-stripe solution of width. Then eigen functions of C (redefined for the 
A^-stripe solution), are close to ^q. Smooth connection of at each boundary x = -^j 
is conditioned. Consideration of Q{x) = Du'^'^ + D^'^'^, which is close to 0^-,^ e^^'^-^S^P, 
gives a rough estimation of the eigen values as A" ~ —4 (^)^c}sPsin^ iw)- '^^^^ indicates 
that a shorter wave (i.e. closer k to y) has larger instability if dsP < 0. A'^ for even 
is identical to that estimated from two stripes. Approximated values are shown in Fig. 01 
The approximated A^ for A^ = 2 is plotted in (a). Eigen values and corresponding eigen 
functions for N = 8, L = 80.0 are shown in (b, c). 

To consider transient processes, an illustrative example can be seen from the stationary 
state of 2A^-stripes with a small perturbation. If Eq. ()14j) is satisfied for wo, the most 
unstable function is \E'^. By the growth of the perturbation along this function, the system 
reaches A^-stripe pattern at last. If this new state becomes unstable, a similar process 
follows until the system reaches a steady state. In the (unstable) stationary state where 
the dynamics become close in their transient, A is small if the distance between adjacent 
stripes is large. The corresponding state lasts for the duration of approximately A~^ and 
therefore each state appears quasi- stationary. Because the distance becomes twice as large 
after each transient, the staying time in the quasi-stationary state also gets longer. We could 
numerically observe these processes from the 8-stripe initial condition. This demonstrates 
the underlying processes of the characteristic transients inherent to mass-conserved systems. 

After the long transient, the system reaches the steady state at which the condition Eq. 
fll4|l is violated. Thus, the characteristic wavelength of the steady state is much longer than 
expected by linear stability analysis at a uniform state in a mass-conserved system. Our 
numerical model showed a one-stripe solution eventually between 2.0 < L < 100.0, while 
£m. = 3.2 (data not shown). 

Notice that Eq. (jSJ, the condition for Turing instability, implies dgP < at the uniform 
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state. Therefore Eq. ()14|) is always satisfied in the early stages of transient, where w{x) 
ranges in the neighborhood of w*. 

One interesting question that arises is the possibility of the system in which the condition 
in Eq. (fT^ is always valid in the transient quasi-steady states except in sole stripe solutions. 
Such systems fall into a one-stripe solution after a long transient regardless of system size. 
An example of such system was supplied by our group in defined by f{u,v) = —a{u + 
v) {{5u + v){u + v) — 13) where 5 = Du/D^. In this specific model, if L is large enough, P is 
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well approximated by P = ^j^^ with 7 = Thus a one-stripe pattern is the 

only stable state in the system. Though rigorous conditions are not described here, many 
mass-conserved systems have such properties. 

In this Letter, we study RDS in which the uniform state is destabilized via a Turing 
mechanism and mass {u + v) is conserved. The analysis presented is useful for stationary 
patterns in any RDS and the conserved quantity does not have to be strictly defined by 
mass. Because the existence of any conserved quantity brings the corresponding 0-eigen 
function, our arguments are applicable. Thus, the dynamics studied here may be observed 
in a wider class of RDS with conserved quantities 

We did not mention the hierarchical structure of quasi-stationary states in the phase- 
space, which is a necessary condition for the sequential transient as discussed in It 
is a global property of the phase space of the systems and difficult to study. Numerical 
simulations suggest it is satisfied in mass-conserved systems. 

Applications of this work are possible to many phenomena, particularly to biological 
systems. Proposed biological models often contain conserved quantities 
the cellular level (~ 10/im), cytosolic proteins diffuse at ~ lO/im^/sec [l7| leading to the 
rough estimation of the time scale of dynamics as A^^ ~ (AP/L^As) ^ ~ LP' /Dy ~ lOsec. 
Typically, it is faster than the synthesis or degradation of molecules and the dynamics is 
expected to occur within the time scale in which mass-conserved modeling is valid. The 
formation of cell polarity based on the above discussions is a potential application 1^ . 
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FIG. 1: Transient dynamics of a mass-conserved system. The model system and the meaning of 
alphabet letters are explained in the text. System size was chosen as L = 50.0 here. Note that time 
scale is represented by log scale. In the system, the most unstable wavelength at homogeneous 
state is £m = 3.2, shown by the line segment in the left bottom. We checked that the system 
eventually falls to a one-stripe pattern for any system size between 2.0 < L < 100.0. 



10 



(a) 



1.0 



0.1 



^ 0.01 



0.001 



0.0001 




(b) 



L=20 





1 t=2500 1 











20 ^0 




1000 



2000 

time 



3000 



4000 



1 \ / \ 
1 \ / I 
1 ^-^ \ 
1 1 
1 1 
/ \ 


' Au 

1 Av 




1 I 1 - 




1 1 ' 




1 1 / V 1 







10 
space X 





I \ / \ du/ds — 
j 1 dv/ds ---- 
1 1 






) \ 
i \ 









space X 


1 



20 



FIG. 2: Two identical steady solutions for ^ = 10.0 are connected and perturbed (±1.0%, keeping 
total mass quantity) at f = 0. (a) |As| grows exponentially with time, indicating that one stripe 
decays while the other grows. u{x,t) (solid) and v{x,t) (broken) at t = 0,2500,3000 and 3500 are 
shown in insets, (b) Aw(x), the difference of w{x,t) between t = 200 and 300 is shown in the top 
panel, while dgWQ is shown in the bottom panel (normalization is applied). 
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FIG. 3: (a) Growth rates A for respective system sizes are evaluated from the observations of 
|As(t)| (x) and from Eq. (fTl^ (O)- Approximated estimation of A, —^dgP, are also plotted 
(+). (b) Eigen values of C for the 8-stripe state are numerically calculated (L = 80.0), and (c) 
corresponding eigen functions. Note that eigen values except k = 4, 8 are degenerated. Sinusoidal 
curves are also shown for guidance. 
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